SUBROUTINE DG_MOT_R
!-----------------------------------------------------------------------
!
! Definition von G:  Bewegungs- bzw. Impulsgleichung
!

      use primvar,  only : DG, X, XZ, XA, MR, MD, MUr, MUp, ME, MH, &
                           DG_i, DG_im1, DG_ip1, DG_im2, DG_ip2
      use physco,   only : z0, z1, z2, z12, pi, sqrt2pi, grav, clight
      use config,   only : np, BCflag
      use global,   only : tst, zz, relax
      use geomvar,  only : V_flux, V_volZ, V_vol, V_volA, S_vol, S_volZ
      use advecvar, only : momr_adv, dmomr_adv_dim1, dmomr_adv_di, dmomr_adv_dip1, dmomr_adv_dip2
      use matvar,   only : Pgas0Z, dPgas0ZdE, dPgas0ZdD, OPArosZ, dOPArosZdE, dOPArosZdD
      use viscvar,  only : divU, divU_flag, alpha_r, q_quad_r, muQ_r, Uterm_r
      use RBvar,    only : M_star
      use zvar,     only : H_pZ, dH_pZdE, dH_pZdR, H_p, dH_pdE, dH_pdR


      implicit none

      integer :: i


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

   ! innere Pseudozellen: i=1 & i=2
      DG(MUr,MUr,DG_i,1)      =  z1
      DG(MUr,MUr,DG_i,2)      =  z1

   ! aeussere Pseudozellen: i=np, i=np-1
      DG(MUr,MUr,DG_i,np)     =  z1
      DG(MUr,MUr,DG_i,np-1)   =  z1
      DG(MUr,MUr,DG_i,np-2)   =  z1


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=3,np-3

DG(MUr,MR,DG_im1,i)=-(pi*X(MD,i)*(X(MR,i-1) + X(MR,i+1))*X(MUr,i)*z12 ) + pi*X(MD,i)*X(MUr,i)*(X(MR,i+1) - &
X(MR,i-1))*z12  - (grav*M_star*pi*XZ(MD,i)*(XZ(MR,i-1) + XZ(MR,i+1))*tst*z12 *zz ) / XZ(MR,i)**2 + &
(grav*M_star*pi*XZ(MD,i)*tst*(XZ(MR,i+1) - XZ(MR,i-1))*z12 *zz ) / XZ(MR,i)**2 - (XZ(MD,i)*tst*z12 *(-(pi*(XZ(MR,i) + &
XZ(MR,i-1))*XZ(MUp,i-1)**2*zz ) + pi*XZ(MUp,i-1)**2*(XZ(MR,i) - XZ(MR,i-1))*zz )) / XZ(MR,i) - (-(pi*X(MR,i-1)*z12 *z2 ) + &
pi*tst*XZ(MUr,i-1)*zz )*momr_adv(i-1) - (XZ(MH,i)*pi*tst*z2 *(-(pi*XZ(MD,i-1)*(XZ(MR,i) + XZ(MR,i-1))*zz *OPArosZ(i-1)) + &
pi*XZ(MD,i-1)*(XZ(MR,i) - XZ(MR,i-1))*zz *OPArosZ(i-1) - (X(MD,i-1)*XZ(MD,i-1)*S_volZ(i-1)*dH_pdR(i-1)*dOPArosZdD(i-1)) / &
(sqrt2pi*H_p(i-1)**2))) / clight + pi*XZ(MR,i)*sqrt2pi*tst*z2 *(-(Pgas0Z(i-1)*dH_pZdR(i-1)) + &
(X(MD,i-1)*H_pZ(i-1)*dH_pdR(i-1)*dPgas0ZdD(i-1)) / (sqrt2pi*H_p(i-1)**2)) - (pi*tst*z12 *(-(muQ_r(i-1)*XZ(MR,i-1)*Uterm_r(i-1)*z2 &
*zz ) - muQ_r(i-1)*(XZ(MR,i)**2 + XZ(MR,i-1)**2)*(z1  - BCflag(i-1))*(-((XZ(MUr,i-1)*zz ) / (XZ(MR,i)**2 + XZ(MR,i-1)**2)) + &
((XZ(MUr,i) - XZ(MUr,i-1))*zz ) / (-(XZ(MR,i-1)) + XZ(MR,i)*(BCflag(i-1) + z1 ))**2 + (XZ(MR,i-1)*(XZ(MR,i)*XZ(MUr,i) + &
XZ(MR,i-1)*XZ(MUr,i-1))*z2 *zz ) / (XZ(MR,i)**2 + XZ(MR,i-1)**2)**2) - (XZ(MR,i)**2 + &
XZ(MR,i-1)**2)*Uterm_r(i-1)*(-(divU_flag(i-1)*divU(i-1)*q_quad_r**2*XZ(MD,i-1)*XZ(MR,i-1)*z2 *zz ) + &
(divU_flag(i-1)*pi*q_quad_r**2*XZ(MD,i-1)*XZ(MR,i-1)**2*XZ(MUr,i-1)*z2 *zz ) / (BCflag(i-1) + S_volZ(i-1)) + &
(divU_flag(i-1)*pi*q_quad_r**2*XZ(MD,i-1)*XZ(MR,i-1)**2*(XZ(MR,i)*XZ(MUr,i) - XZ(MR,i-1)*XZ(MUr,i-1))*z2 *(-(pi*(XZ(MR,i) + &
XZ(MR,i-1))*zz ) + pi*(XZ(MR,i) - XZ(MR,i-1))*zz )) / (BCflag(i-1) + S_volZ(i-1))**2 + &
(alpha_r(i-1)*XZ(MD,i-1)*Sqrt(Pgas0Z(i-1))*dH_pZdR(i-1)) / Sqrt(XZ(MD,i-1) / (sqrt2pi*H_pZ(i-1))) + &
(alpha_r(i-1)*XZ(MD,i-1)**2*z12 *Sqrt(Pgas0Z(i-1))*dH_pZdR(i-1)) / (sqrt2pi*(XZ(MD,i-1) / (sqrt2pi*H_pZ(i-1)))**1.5*H_pZ(i-1)) - &
(alpha_r(i-1)*X(MD,i-1)*XZ(MD,i-1)*z12 *H_pZ(i-1)*dH_pdR(i-1)*dPgas0ZdD(i-1)) / (sqrt2pi*H_p(i-1)**2*Sqrt(XZ(MD,i-1) / &
(sqrt2pi*H_pZ(i-1)))*Sqrt(Pgas0Z(i-1)))))) / XZ(MR,i)
DG(MUr,MR,DG_i,i)=(XZ(MD,i)*tst*(S_volZ(i)*XZ(MUp,i)**2 + S_volZ(i-1)*XZ(MUp,i-1)**2)*z12 *zz ) / XZ(MR,i)**2 + &
(pi*tst*(muQ_r(i)*(XZ(MR,i)**2 + XZ(MR,i+1)**2)*Uterm_r(i) - muQ_r(i-1)*(XZ(MR,i)**2 + XZ(MR,i-1)**2)*Uterm_r(i-1))*z12 *zz ) / &
XZ(MR,i)**2 - (grav*M_star*XZ(MD,i)*tst*V_volZ(i)*z2 *zz ) / XZ(MR,i)**3 - (XZ(MD,i)*tst*z12 *(pi*(XZ(MR,i) + &
XZ(MR,i-1))*XZ(MUp,i-1)**2*zz  - pi*(XZ(MR,i) + XZ(MR,i+1))*XZ(MUp,i)**2*zz  + pi*XZ(MUp,i)**2*(XZ(MR,i+1) - XZ(MR,i))*zz  + &
pi*XZ(MUp,i-1)**2*(XZ(MR,i) - XZ(MR,i-1))*zz )) / XZ(MR,i) + (-(pi*X(MR,i)*z12 *z2 ) + pi*tst*XZ(MUr,i)*zz )*momr_adv(i) - &
(-(pi*X(MR,i)*z12 *z2 ) + pi*tst*XZ(MUr,i)*zz )*momr_adv(i-1) + pi*sqrt2pi*tst*z2 *zz *(H_pZ(i)*Pgas0Z(i) - &
H_pZ(i-1)*Pgas0Z(i-1)) - (XZ(MH,i)*pi*tst*z2 *(-(pi*XZ(MD,i)*(XZ(MR,i) + XZ(MR,i+1))*zz *OPArosZ(i)) + pi*XZ(MD,i)*(XZ(MR,i+1) - &
XZ(MR,i))*zz *OPArosZ(i) + pi*XZ(MD,i-1)*(XZ(MR,i) + XZ(MR,i-1))*zz *OPArosZ(i-1) + pi*XZ(MD,i-1)*(XZ(MR,i) - XZ(MR,i-1))*zz &
*OPArosZ(i-1) - (X(MD,i)*XZ(MD,i)*S_volZ(i)*dH_pdR(i)*dOPArosZdD(i)) / (sqrt2pi*H_p(i)**2))) / clight + &
pi*XZ(MR,i)*sqrt2pi*tst*z2 *(Pgas0Z(i)*dH_pZdR(i) - (X(MD,i)*H_pZ(i)*dH_pdR(i)*dPgas0ZdD(i)) / (sqrt2pi*H_p(i)**2)) - (pi*tst*z12 &
*(muQ_r(i)*XZ(MR,i)*Uterm_r(i)*z2 *zz  - muQ_r(i-1)*XZ(MR,i)*Uterm_r(i-1)*z2 *zz  - muQ_r(i-1)*(XZ(MR,i)**2 + XZ(MR,i-1)**2)*(z1  &
- BCflag(i-1))*(-((XZ(MUr,i)*zz ) / (XZ(MR,i)**2 + XZ(MR,i-1)**2)) - ((BCflag(i-1) + z1 )*(XZ(MUr,i) - XZ(MUr,i-1))*zz ) / &
(-(XZ(MR,i-1)) + XZ(MR,i)*(BCflag(i-1) + z1 ))**2 + (XZ(MR,i)*(XZ(MR,i)*XZ(MUr,i) + XZ(MR,i-1)*XZ(MUr,i-1))*z2 *zz ) / &
(XZ(MR,i)**2 + XZ(MR,i-1)**2)**2) + muQ_r(i)*(XZ(MR,i)**2 + XZ(MR,i+1)**2)*(z1  - BCflag(i))*(-((XZ(MUr,i)*zz ) / (XZ(MR,i)**2 + &
XZ(MR,i+1)**2)) + ((XZ(MUr,i+1) - XZ(MUr,i))*zz ) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 ))**2 + &
(XZ(MR,i)*(XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))*z2 *zz ) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)**2) - (XZ(MR,i)**2 + &
XZ(MR,i-1)**2)*Uterm_r(i-1)*(-((divU_flag(i-1)*pi*q_quad_r**2*XZ(MD,i-1)*XZ(MR,i-1)**2*XZ(MUr,i)*z2 *zz ) / (BCflag(i-1) + &
S_volZ(i-1))) + (divU_flag(i-1)*pi*q_quad_r**2*XZ(MD,i-1)*XZ(MR,i-1)**2*(XZ(MR,i)*XZ(MUr,i) - XZ(MR,i-1)*XZ(MUr,i-1))*z2 &
*(pi*(XZ(MR,i) + XZ(MR,i-1))*zz  + pi*(XZ(MR,i) - XZ(MR,i-1))*zz )) / (BCflag(i-1) + S_volZ(i-1))**2) + (XZ(MR,i)**2 + &
XZ(MR,i+1)**2)*Uterm_r(i)*(-(divU_flag(i)*divU(i)*q_quad_r**2*XZ(MD,i)*XZ(MR,i)*z2 *zz ) + &
(divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MUr,i)*z2 *zz ) / (BCflag(i) + S_volZ(i)) + &
(divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 *(-(pi*(XZ(MR,i) + &
XZ(MR,i+1))*zz ) + pi*(XZ(MR,i+1) - XZ(MR,i))*zz )) / (BCflag(i) + S_volZ(i))**2 + &
(alpha_r(i)*XZ(MD,i)*Sqrt(Pgas0Z(i))*dH_pZdR(i)) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + (alpha_r(i)*XZ(MD,i)**2*z12 &
*Sqrt(Pgas0Z(i))*dH_pZdR(i)) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5*H_pZ(i)) - (alpha_r(i)*X(MD,i)*XZ(MD,i)*z12 &
*H_pZ(i)*dH_pdR(i)*dPgas0ZdD(i)) / (sqrt2pi*H_p(i)**2*Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i)))))) / XZ(MR,i)
DG(MUr,MR,DG_ip1,i)=pi*X(MD,i)*(X(MR,i-1) + X(MR,i+1))*X(MUr,i)*z12  + pi*X(MD,i)*X(MUr,i)*(X(MR,i+1) - &
X(MR,i-1))*z12  + (grav*M_star*pi*XZ(MD,i)*(XZ(MR,i-1) + XZ(MR,i+1))*tst*z12 *zz ) / XZ(MR,i)**2 + &
(grav*M_star*pi*XZ(MD,i)*tst*(XZ(MR,i+1) - XZ(MR,i-1))*z12 *zz ) / XZ(MR,i)**2 - (XZ(MD,i)*tst*z12 *(pi*(XZ(MR,i) + &
XZ(MR,i+1))*XZ(MUp,i)**2*zz  + pi*XZ(MUp,i)**2*(XZ(MR,i+1) - XZ(MR,i))*zz )) / XZ(MR,i) - (pi*tst*z12 &
*(muQ_r(i)*XZ(MR,i+1)*Uterm_r(i)*z2 *zz  + muQ_r(i)*(XZ(MR,i)**2 + XZ(MR,i+1)**2)*(z1  - BCflag(i))*(-((XZ(MUr,i+1)*zz ) / &
(XZ(MR,i)**2 + XZ(MR,i+1)**2)) - ((BCflag(i) + z1 )*(XZ(MUr,i+1) - XZ(MUr,i))*zz ) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 &
))**2 + (XZ(MR,i+1)*(XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))*z2 *zz ) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)**2) + (XZ(MR,i)**2 + &
XZ(MR,i+1)**2)*Uterm_r(i)*(-((divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MUr,i+1)*z2 *zz ) / (BCflag(i) + S_volZ(i))) + &
(divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 *(pi*(XZ(MR,i) + &
XZ(MR,i+1))*zz  + pi*(XZ(MR,i+1) - XZ(MR,i))*zz )) / (BCflag(i) + S_volZ(i))**2))) / XZ(MR,i) + (-(pi*X(MR,i+1)*z12 *z2 ) + &
pi*tst*XZ(MUr,i+1)*zz )*momr_adv(i) - (XZ(MH,i)*pi*tst*z2 *(pi*XZ(MD,i)*(XZ(MR,i) + XZ(MR,i+1))*zz *OPArosZ(i) + &
pi*XZ(MD,i)*(XZ(MR,i+1) - XZ(MR,i))*zz *OPArosZ(i))) / clight
DG(MUr,MD,DG_im2,i)=-(XZ(MUr,i-2)*V_flux(i-1)*zz *dmomr_adv_dim1(i-1))
DG(MUr,MD,DG_im1,i)=-((XZ(MH,i)*pi*tst*z2 *(S_volZ(i-1)*zz *OPArosZ(i-1) + &
(XZ(MD,i-1)*S_volZ(i-1)*dOPArosZdD(i-1)) / (sqrt2pi*H_p(i-1)))) / clight) - (pi*XZ(MR,i)*tst*z2 *H_pZ(i-1)*dPgas0ZdD(i-1)) / &
H_p(i-1) + (pi*(XZ(MR,i)**2 + XZ(MR,i-1)**2)*tst*Uterm_r(i-1)*z12 *(-(divU_flag(i-1)*divU(i-1)*q_quad_r**2*XZ(MR,i-1)**2*zz ) - &
(alpha_r(i-1)*XZ(MD,i-1)*z12 *zz *Sqrt(Pgas0Z(i-1))) / (sqrt2pi*(XZ(MD,i-1) / (sqrt2pi*H_pZ(i-1)))**1.5) + (alpha_r(i-1)*zz &
*H_pZ(i-1)*Sqrt(Pgas0Z(i-1))) / Sqrt(XZ(MD,i-1) / (sqrt2pi*H_pZ(i-1))) + (alpha_r(i-1)*XZ(MD,i-1)*z12 *H_pZ(i-1)*dPgas0ZdD(i-1)) &
/ (sqrt2pi*H_p(i-1)*Sqrt(XZ(MD,i-1) / (sqrt2pi*H_pZ(i-1)))*Sqrt(Pgas0Z(i-1))))) / XZ(MR,i) - XZ(MUr,i-1)*V_flux(i-1)*zz &
*dmomr_adv_di(i-1) + XZ(MUr,i-1)*V_flux(i)*zz *dmomr_adv_dim1(i)
DG(MUr,MD,DG_i,i)=X(MUr,i)*V_vol(i) + (grav*M_star*tst*V_volZ(i)*zz ) / XZ(MR,i)**2 - &
(tst*(S_volZ(i)*XZ(MUp,i)**2 + S_volZ(i-1)*XZ(MUp,i-1)**2)*z12 *zz ) / XZ(MR,i) - (XZ(MH,i)*pi*tst*z2 *(S_volZ(i)*zz *OPArosZ(i) &
+ (XZ(MD,i)*S_volZ(i)*dOPArosZdD(i)) / (sqrt2pi*H_p(i)))) / clight + (pi*XZ(MR,i)*tst*z2 *H_pZ(i)*dPgas0ZdD(i)) / H_p(i) - &
(pi*(XZ(MR,i)**2 + XZ(MR,i+1)**2)*tst*Uterm_r(i)*z12 *(-(divU_flag(i)*divU(i)*q_quad_r**2*XZ(MR,i)**2*zz ) - &
(alpha_r(i)*XZ(MD,i)*z12 *zz *Sqrt(Pgas0Z(i))) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5) + (alpha_r(i)*zz &
*H_pZ(i)*Sqrt(Pgas0Z(i))) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + (alpha_r(i)*XZ(MD,i)*z12 *H_pZ(i)*dPgas0ZdD(i)) / &
(sqrt2pi*H_p(i)*Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i))))) / XZ(MR,i) - XZ(MUr,i)*V_flux(i-1)*zz *dmomr_adv_dip1(i-1) &
+ XZ(MUr,i)*V_flux(i)*zz *dmomr_adv_di(i)
DG(MUr,MD,DG_ip1,i)=-(XZ(MUr,i+1)*V_flux(i-1)*zz *dmomr_adv_dip2(i-1)) + XZ(MUr,i+1)*V_flux(i)*zz &
*dmomr_adv_dip1(i)
DG(MUr,MD,DG_ip2,i)=XZ(MUr,i+2)*V_flux(i)*zz *dmomr_adv_dip2(i)
DG(MUr,ME,DG_im1,i)=-((XZ(MH,i)*pi**2*XZ(MD,i-1)*(XZ(MR,i) + XZ(MR,i-1))*tst*(XZ(MR,i) - XZ(MR,i-1))*z2 &
*(-((X(MD,i-1)*dOPArosZdD(i-1)*dH_pdE(i-1)) / (sqrt2pi*H_p(i-1)**2)) + dOPArosZdE(i-1))) / clight) + pi*XZ(MR,i)*sqrt2pi*tst*z2 &
*(-(Pgas0Z(i-1)*dH_pZdE(i-1)) - H_pZ(i-1)*(-((X(MD,i-1)*dPgas0ZdD(i-1)*dH_pdE(i-1)) / (sqrt2pi*H_p(i-1)**2)) + dPgas0ZdE(i-1))) + &
(pi*(XZ(MR,i)**2 + XZ(MR,i-1)**2)*tst*Uterm_r(i-1)*z12 *((alpha_r(i-1)*XZ(MD,i-1)*Sqrt(Pgas0Z(i-1))*dH_pZdE(i-1)) / &
Sqrt(XZ(MD,i-1) / (sqrt2pi*H_pZ(i-1))) + (alpha_r(i-1)*XZ(MD,i-1)**2*z12 *Sqrt(Pgas0Z(i-1))*dH_pZdE(i-1)) / (sqrt2pi*(XZ(MD,i-1) &
/ (sqrt2pi*H_pZ(i-1)))**1.5*H_pZ(i-1)) + (alpha_r(i-1)*XZ(MD,i-1)*z12 *H_pZ(i-1)*(-((X(MD,i-1)*dPgas0ZdD(i-1)*dH_pdE(i-1)) / &
(sqrt2pi*H_p(i-1)**2)) + dPgas0ZdE(i-1))) / (Sqrt(XZ(MD,i-1) / (sqrt2pi*H_pZ(i-1)))*Sqrt(Pgas0Z(i-1))))) / XZ(MR,i)
DG(MUr,ME,DG_i,i)=-((XZ(MH,i)*pi**2*XZ(MD,i)*(XZ(MR,i) + XZ(MR,i+1))*tst*(XZ(MR,i+1) - XZ(MR,i))*z2 &
*(-((X(MD,i)*dOPArosZdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dOPArosZdE(i))) / clight) + pi*XZ(MR,i)*sqrt2pi*tst*z2 &
*(Pgas0Z(i)*dH_pZdE(i) + H_pZ(i)*(-((X(MD,i)*dPgas0ZdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dPgas0ZdE(i))) - (pi*(XZ(MR,i)**2 + &
XZ(MR,i+1)**2)*tst*Uterm_r(i)*z12 *((alpha_r(i)*XZ(MD,i)*Sqrt(Pgas0Z(i))*dH_pZdE(i)) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + &
(alpha_r(i)*XZ(MD,i)**2*z12 *Sqrt(Pgas0Z(i))*dH_pZdE(i)) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5*H_pZ(i)) + &
(alpha_r(i)*XZ(MD,i)*z12 *H_pZ(i)*(-((X(MD,i)*dPgas0ZdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dPgas0ZdE(i))) / (Sqrt(XZ(MD,i) / &
(sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i))))) / XZ(MR,i)
DG(MUr,MUr,DG_im2,i)=-(XZ(MD,i-2)*V_flux(i-1)*zz *dmomr_adv_dim1(i-1))
DG(MUr,MUr,DG_im1,i)=-((pi*tst*z12 *(-((divU_flag(i-1)*pi*q_quad_r**2*XZ(MD,i-1)*XZ(MR,i-1)**3*(XZ(MR,i)**2 + &
XZ(MR,i-1)**2)*Uterm_r(i-1)*z2 *zz ) / (BCflag(i-1) + S_volZ(i-1))) - muQ_r(i-1)*(XZ(MR,i)**2 + XZ(MR,i-1)**2)*(z1  - &
BCflag(i-1))*(-((XZ(MR,i-1)*zz ) / (XZ(MR,i)**2 + XZ(MR,i-1)**2)) - (zz ) / (-(XZ(MR,i-1)) + XZ(MR,i)*(BCflag(i-1) + z1 ))))) / &
XZ(MR,i)) - pi*XZ(MR,i-1)*tst*zz *momr_adv(i-1) - XZ(MD,i-1)*V_flux(i-1)*zz *dmomr_adv_di(i-1) + XZ(MD,i-1)*V_flux(i)*zz &
*dmomr_adv_dim1(i)
DG(MUr,MUr,DG_i,i)=X(MD,i)*V_vol(i) - (pi*tst*z12 &
*((divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**3*(XZ(MR,i)**2 + XZ(MR,i+1)**2)*Uterm_r(i)*z2 *zz ) / (BCflag(i) + S_volZ(i)) + &
(divU_flag(i-1)*pi*q_quad_r**2*XZ(MD,i-1)*XZ(MR,i)*XZ(MR,i-1)**2*(XZ(MR,i)**2 + XZ(MR,i-1)**2)*Uterm_r(i-1)*z2 *zz ) / &
(BCflag(i-1) + S_volZ(i-1)) + muQ_r(i)*(XZ(MR,i)**2 + XZ(MR,i+1)**2)*(z1  - BCflag(i))*(-((XZ(MR,i)*zz ) / (XZ(MR,i)**2 + &
XZ(MR,i+1)**2)) - (zz ) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 ))) - muQ_r(i-1)*(XZ(MR,i)**2 + XZ(MR,i-1)**2)*(z1  - &
BCflag(i-1))*(-((XZ(MR,i)*zz ) / (XZ(MR,i)**2 + XZ(MR,i-1)**2)) + zz  / (-(XZ(MR,i-1)) + XZ(MR,i)*(BCflag(i-1) + z1 ))))) / &
XZ(MR,i) + pi*XZ(MR,i)*tst*zz *momr_adv(i) - pi*XZ(MR,i)*tst*zz *momr_adv(i-1) - XZ(MD,i)*V_flux(i-1)*zz *dmomr_adv_dip1(i-1) + &
XZ(MD,i)*V_flux(i)*zz *dmomr_adv_di(i)
DG(MUr,MUr,DG_ip1,i)=-((pi*tst*z12 *(-((divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MR,i+1)*(XZ(MR,i)**2 &
+ XZ(MR,i+1)**2)*Uterm_r(i)*z2 *zz ) / (BCflag(i) + S_volZ(i))) + muQ_r(i)*(XZ(MR,i)**2 + XZ(MR,i+1)**2)*(z1  - &
BCflag(i))*(-((XZ(MR,i+1)*zz ) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + zz  / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 ))))) / &
XZ(MR,i)) + pi*XZ(MR,i+1)*tst*zz *momr_adv(i) - XZ(MD,i+1)*V_flux(i-1)*zz *dmomr_adv_dip2(i-1) + XZ(MD,i+1)*V_flux(i)*zz &
*dmomr_adv_dip1(i)
DG(MUr,MUr,DG_ip2,i)=XZ(MD,i+2)*V_flux(i)*zz *dmomr_adv_dip2(i)
DG(MUr,MUp,DG_im1,i)=-((XZ(MD,i)*S_volZ(i-1)*tst*XZ(MUp,i-1)*z12 *z2 *zz ) / XZ(MR,i))
DG(MUr,MUp,DG_i,i)=-((XZ(MD,i)*S_volZ(i)*tst*XZ(MUp,i)*z12 *z2 *zz ) / XZ(MR,i))
DG(MUr,MH,DG_i,i)=-((pi*tst*z2 *zz *(XZ(MD,i)*S_volZ(i)*OPArosZ(i) + XZ(MD,i-1)*S_volZ(i-1)*OPArosZ(i-1))) / &
clight)

      end do


!-----------------------------------------------------------------------
!    Geschlossene Raender fuers Relaxieren
!-----------------------------------------------------------------------

! Ueberschreibt einfach Randwerte aus obiger Schleife
      DG(MUr,:,:,np-2)        =  z0                                !debugging
      DG(MUr,MUr,DG_i,np-2)   =  z1                                !debugging
      if (relax) then                                              !debugging
      ! innere Randbedingung: i=3
         DG(MUr,:,:,3)           =  z0
         DG(MUr,MUr,DG_i,3)      =  z1
      ! aeussere Randbedingung: i=np-2
         DG(MUr,:,:,np-2)        =  z0
         DG(MUr,MUr,DG_i,np-2)   =  z1
      end if


END SUBROUTINE DG_MOT_R
